Introduction

In this example, we deal with effects of attendance on final exam performance

A regression model is used to explain the outcome on a final exam (stndfnl, standardized observations) as a function of the following regressors:

  • percentage of classes attended,
  • prior college grade point average,
  • ACT score (American College Testing).

Data:

attend <- load("attend.RData")
print(desc)
##    variable                        label
## 1    attend   classes attended out of 32
## 2   termGPA                 GPA for term
## 3    priGPA cumulative GPA prior to term
## 4       ACT                    ACT score
## 5     final             final exam score
## 6   atndrte     percent classes attended
## 7     hwrte   percent homework turned in
## 8     frosh               =1 if freshman
## 9      soph              =1 if sophomore
## 10   missed     number of classes missed
## 11  stndfnl            (final - mean)/sd

Base model (no interactions)

\[\textit{stndfnl}_{\,i} = \beta_0 + \beta_1 \textit{atndrte}_i + \beta_2 \textit{priGPA}_i + \beta_3 \textit{ACT}_i + u_i\]

lm1<-lm(stndfnl ~ atndrte + priGPA + ACT, data=data)
summary(lm1)
## 
## Call:
## lm(formula = stndfnl ~ atndrte + priGPA + ACT, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2339 -0.5528 -0.0329  0.5884  2.3303 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.343655   0.299099 -11.179  < 2e-16 ***
## atndrte      0.005334   0.002369   2.252   0.0247 *  
## priGPA       0.402373   0.078280   5.140 3.60e-07 ***
## ACT          0.084257   0.011182   7.535 1.57e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.8862 on 676 degrees of freedom
## Multiple R-squared:  0.2013, Adjusted R-squared:  0.1978 
## F-statistic: 56.79 on 3 and 676 DF,  p-value: < 2.2e-16

The above standard errors come from the OLS estimation of variance-covariance matrix of estimated parameters: \[ \textit{var}(\hat{\mathbf{\beta}})=\hat{\sigma}^2(\mathbf{X}^{\prime}\mathbf{X})^{-1}. \]
Different types of robust standard errors may be calculated, White’s HC0 are given by: \[ \textit{var}(\hat{\mathbf{\beta}})= \mathbf{\hat{\,\Sigma \,}} = [(\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{X}^{\prime}] \mathbf{\hat{\,\Omega \,}} [(\mathbf{X}^{\prime}\mathbf{X})^{-1}\mathbf{X}^{\prime}]^{\prime} \] where \(\mathbf{\hat{\Omega}}\) is a diagonal matrix with \((i,i)\) elements equal to \(\hat{u}_i^2\).

coeftest(lm1, vcov = vcovHC(lm1, "HC0")) # White's HC0 robust errors
## 
## t test of coefficients:
## 
##               Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept) -3.3436549  0.2922276 -11.4420 < 2.2e-16 ***
## atndrte      0.0053337  0.0023459   2.2736    0.0233 *  
## priGPA       0.4023727  0.0831316   4.8402 1.609e-06 ***
## ACT          0.0842571  0.0110510   7.6244 8.331e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Here, given the linear nature of the model, the estimated \(\hat{\beta}_1\) to \(\hat{\beta}_3\) coefficients are:

  • the ceteris paribus marginal effects. They are valid regardless of magnitudes of their corresponding regressors or magnitudes of other regressors.
    • Can be viewed as Average partial effects, APEs (average marginal effects)
    • Can be viewed as Partial effects at the average, PEAs (marginal effects at the average)

Standard errors \(\textit{s.e.}(\hat{\beta}_j)\) are valid, regardless of observed regressor values (magnitudes).

Model with interaction term: atndrte * priGPA

This LRM may be outlined as:

\[\textit{stndfnl}_{\,i} = \beta_0 + \beta_1 \textit{atndrte}_i + \beta_2 \textit{priGPA}_i + \beta_3 \textit{ACT}_i + \beta_4 (\textit{atndrte}_i \times \textit{priGPA}_i ) + u_i\]

lm2<-lm(stndfnl ~ atndrte + priGPA + ACT + atndrte:priGPA , data=data)
summary(lm2)
## 
## Call:
## lm(formula = stndfnl ~ atndrte + priGPA + ACT + atndrte:priGPA, 
##     data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2071 -0.5380 -0.0297  0.5852  2.3765 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    -1.135889   0.793175  -1.432  0.15258    
## atndrte        -0.020893   0.009047  -2.309  0.02122 *  
## priGPA         -0.554498   0.328065  -1.690  0.09145 .  
## ACT             0.081698   0.011149   7.328  6.7e-13 ***
## atndrte:priGPA  0.011462   0.003818   3.002  0.00278 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.881 on 675 degrees of freedom
## Multiple R-squared:  0.2118, Adjusted R-squared:  0.2072 
## F-statistic: 45.35 on 4 and 675 DF,  p-value: < 2.2e-16

Now, the partial effects of atndrte, priGPA and ACT are as follows:

  • \(\frac{\partial \textit{stndfnl}}{\partial \textit{atndrte}_i} = \beta_1 + \beta_4 \textit{priGPA}_i\)

  • \(\frac{\partial \textit{stndfnl}}{\partial \textit{priGPA}_i} = \beta_2 + \beta_4 \textit{atndrte}_i\)

  • \(\frac{\partial \textit{stndfnl}}{\partial \textit{ACT}} = \beta_3\)

Partial effects:

  • \(\textit{ACT}\): Interpretation of \(\beta_3\) does not change from previous example (no interaction).

  • \(\textit{atndrte}\): The partial effect on the dependent variable with respect to \(\textit{atndrte}\) depends on the magnitude of \(\textit{priGPA}\).

    • \(\beta_1\) and \(\textit{s.e.}(\hat{\beta}_1)\) are valid (directly applicable) only if \(\textit{priGPA} =0\). Generally, this is not very useful.
    • For any \(\textit{priGPA} \neq 0\), the partial effect of \(\textit{atndrte}\) is easy to calculate. However, standard error of the \(\textit{atndrte}\) partial effect is not directly available here.
    • Reparametrization and Delta method may be used to get the s.e. (Bootstrap as well).

  • \(\textit{priGPA}\): same situation as above, partial effects depend on the magnitude of the interacting element.

Marginal effects - Reparametrization

Say, we are interested in partial effects and standard errors for a representative student with average atndrte and priGPA values. The effects can be calculated easily:

m1 <- mean(data$atndrte)
m2 <- mean(data$priGPA)
m1 # mean(data$atndrte)
## [1] 81.70956
m2 # mean(data$priGPA)
## [1] 2.586775

  • \(\frac{\Delta \textit{stndfnl}}{\Delta \textit{atndrte}_i} \doteq -0.021 + 0.011 \, \overline{\textit{priGPA}} = -0.021 + 0.011 \times 2.587 = 0.009\)

  • \(\frac{\Delta \textit{stndfnl}}{\Delta \textit{priGPA}_i} \doteq - 0.555 + 0.011 \, \overline{\textit{atndrte}} = - 0.555 + 0.011 \times 81.710 = 0.382\)

  • Standard errors to these effects are available through reparametrization of the interaction term: \(([\textit{atndrte}_i - \overline{\textit{atndrte}} ] \times [\textit{priGPA}_i - \overline{\textit{priGPA}} ] )\)

  • In the reparametrized model, \(\beta_1\) is the effect of unit change in \(\textit{atndrte}\) while \(\textit{priGPA} = \overline{\textit{priGPA}}\)

  • \(\beta_2\) can be interpreted by analogy (in the reparametrized model).

lm3<-lm(stndfnl ~ atndrte + priGPA + ACT + I((atndrte-m1)*(priGPA-m2)) , data=data)
summary(lm3)
## 
## Call:
## lm(formula = stndfnl ~ atndrte + priGPA + ACT + I((atndrte - 
##     m1) * (priGPA - m2)), data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.2071 -0.5380 -0.0297  0.5852  2.3765 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -3.558479   0.305829 -11.636  < 2e-16 ***
## atndrte                            0.008756   0.002616   3.347 0.000863 ***
## priGPA                             0.382031   0.078115   4.891 1.26e-06 ***
## ACT                                0.081698   0.011149   7.328 6.70e-13 ***
## I((atndrte - m1) * (priGPA - m2))  0.011462   0.003818   3.002 0.002778 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.881 on 675 degrees of freedom
## Multiple R-squared:  0.2118, Adjusted R-squared:  0.2072 
## F-statistic: 45.35 on 4 and 675 DF,  p-value: < 2.2e-16

Marginal effects using {margins} package

Again, we are interested in partial effects and standard errors for a representative student with average atndrte and priGPA values. The effects can be calculated easily through {margins} package

lm2<-lm(stndfnl ~ atndrte + priGPA + ACT + atndrte:priGPA , data=data) # repeated
margins(lm2, at = list(atndrte = m1, priGPA = m2)) # see help for {margins} package
## Average marginal effects at specified values
## lm(formula = stndfnl ~ atndrte + priGPA + ACT + atndrte:priGPA,     data = data)
##  at(atndrte) at(priGPA)  atndrte priGPA    ACT
##        81.71      2.587 0.008756  0.382 0.0817

Standard errors are not yet implemented in {margins} for the at= argument.

Delta method - marginal effects & standard errors.

Delta method intuition

Consider a vector of estimated parameters \(\hat{\mathbf{\beta}}\) and some function \(\hat{g} = f(\hat{\mathbf{\beta}})\). Both \(\hat{\mathbf{\beta}}\) and \(\hat{g}\) are random variables and the variance of \(\hat{g}\) is a function of variance of \(\hat{\mathbf{\beta}}\).

  1. If \(g\) is a linear function of \(\mathbf{\beta}\), say \(g = \mathbf{\omega}^{\prime}\mathbf{\beta}\) then: \[ \text{var} (\hat{g}) = \mathbf{\omega}^{\prime} \text{var}(\hat{\mathbf{\beta}}) \mathbf{\omega}\] (recall how variance of fitted value \(\hat{y}_i\) is calculated).

  2. If \(g\) is non-linear, variance of \(\hat{g}\) can be approximated through a first-order Taylor expansion.

Specifically, if \(g\) is a function of parameter \(\beta\) and \(\hat{\beta}\) is a consistent, normally distributed estimator for that parameter: \[ g(\hat{\beta}) \approx g(\beta) + \nabla g(\beta)^\prime (\hat{\beta} - \beta) \] Since \(\beta\) is a constant, and \(\hat{\beta}\) is a consistent estimator for \(\beta\), we can then say: \[ \sqrt{n}\left(g(\hat{\beta})-g(\beta)\right)\,\xrightarrow{D}\,N\left(0, \nabla g(\beta)^\prime \cdot \mathbf{\Sigma}_{\hat{\beta}} \cdot \nabla g(\beta)\right) \]

where \(\hat{\beta}\) and \(\Sigma_{\hat{\beta}}\) come from the OLS estimate and \(g(\mathbf{\beta})\) is the marginal effect.

Emprical example 1, using attend data

In this exmple, we consider a simple Delta method example, as we only deal with linear functions of \(\beta\).

For our linear model,

\[\textit{stndfnl}_{\,i} = \beta_0 + \beta_1 \textit{atndrte}_i + \beta_2 \textit{priGPA}_i + \beta_3 \textit{ACT}_i + \beta_4 (\textit{atndrte}_i \times \textit{priGPA}_i ) + u_i,\]

Let’s consider a unit change in atndrte (and fixing \(\textit{priGPA}=\overline{\textit{priGPA}}\)):

  • Marginal effect can be expressed as a linear function (linear in parameters): \[ g(\mathbf{\beta}) = \beta_1 + \beta_4 \overline{\textit{priGPA}} \]

  • If you take the gradient of this function (remember, a function of relevant coefficients not of regressors), it would be: \[ \nabla g(\mathbf{\beta})^\prime = [1,\,\, \overline{\textit{priGPA}}]^{\,\prime} \,, \]

  • The matrix \(\Sigma_{\hat{\beta}}\) comes from the OLS estimation of our main regression:

\[ \Sigma_{\hat{\beta}} = \left[ \begin{matrix} \textit{var}(\hat{\beta}_1) & \textit{cov}(\hat{\beta}_1,\hat{\beta}_4) \\ \textit{cov}(\hat{\beta}_1,\hat{\beta}_4) & \textit{var}(\hat{\beta}_4) \end{matrix}\right] \]

  • Combining the above, the variance of \(\hat{g} \equiv g(\hat{\beta})\) can be calculated using \(\textit{var} (g) = \nabla g(\beta)^\prime \cdot \mathbf{\Sigma}_{\hat{\beta}} \cdot \nabla g(\beta)\), i.e. as: \[ \textit{var}(\hat{g}) = [1,\,\, \overline{\textit{priGPA}}] \left[ \begin{matrix} \textit{var}(\hat{\beta}_1) & \textit{cov}(\hat{\beta}_1,\hat{\beta}_4) \\ \textit{cov}(\hat{\beta}_1,\hat{\beta}_4) & \textit{var}(\hat{\beta}_4) \end{matrix}\right] \left[ \begin{matrix} 1 \\ \overline{\textit{priGPA}} \end{matrix}\right] \]

Marginal effect & statistical inference of a unit change in atndrte (fixing priGPA at mean) can be calculated as follows:

V <- vcov(lm2)
V # Full VCV matrix for beta coeffs of our LM
##                 (Intercept)       atndrte        priGPA           ACT
## (Intercept)     0.629126678 -6.842976e-03 -0.2330259251 -3.132360e-03
## atndrte        -0.006842976  8.184704e-05  0.0026881773  1.696422e-05
## priGPA         -0.233025925  2.688177e-03  0.1076267846 -1.357526e-04
## ACT            -0.003132360  1.696422e-05 -0.0001357526  1.243008e-04
## atndrte:priGPA  0.002807137 -3.334626e-05 -0.0012166451 -3.253989e-06
##                atndrte:priGPA
## (Intercept)      2.807137e-03
## atndrte         -3.334626e-05
## priGPA          -1.216645e-03
## ACT             -3.253989e-06
## atndrte:priGPA   1.457334e-05
V1 <- V[c(2,5),c(2,5)]
V1 # VCV extracted for beta1 and beta4
##                      atndrte atndrte:priGPA
## atndrte         8.184704e-05  -3.334626e-05
## atndrte:priGPA -3.334626e-05   1.457334e-05
g1 <- c(1, m2) # gradient
g1 # show gradient (as row vector)
## [1] 1.000000 2.586775
VC1 <- t(g1) %*% V1 %*% g1
sqrt(VC1) # s.e. of the marginal effect through delta metod. Compare to "Reparametrization"
##             [,1]
## [1,] 0.002616216

Emprical example 2, using attend data

For a unit change in priGPA (and fixing \(\textit{atndrte}=\overline{\textit{atndrte}}\)):

V <- vcov(lm2) # Full VCV, (repeated for clarity)
V2 <- V[c(3,5),c(3,5)]  # (Sigma matrix for beta2 and beta4)
g2 <- c(1, m1) # gradient
VC2 <- t(g2) %*% V2 %*% g2
sqrt(VC2) # compare to "Reparametrization"
##            [,1]
## [1,] 0.07811469

Emprical example 3, using attend data

The previous two steps can be combined as follows:

Let’s define a matrix \(\mathbf{G}=\{\partial g_i(\mathbf{\beta}) / \partial(\beta_j)\}\)

In our case (for marginal effects of the two interacting regressors):

\[ \mathbf{G}= \left[ \begin{matrix} 1 & 0 & \overline{\textit{priGPA}} \\ 0 & 1 & \overline{\textit{atndrte}} \end{matrix}\right] \]

Now, the variance of marginal effects can be estimated as \[ \textit{var}[g(\mathbf{\hat{\beta}})] = \mathbf{G} \,\,\mathbf{\Sigma}_{\hat{\beta}} \, \mathbf{G}^{\prime} \,, \]

V <- vcov(lm2) # repeated
V3 <- V[c(2,3,5),c(2,3,5)]
G <- matrix(c(1,0,m2,0,1,m1), nrow=2,byrow = T)
VCV <- G %*% V3 %*% t(G)
sqrt(diag(VCV)) # compare to "Reparametrization"
## [1] 0.002616216 0.078114686

Delta method using the {car} and {RcmdrMisc} packages

(Please note that DeltaMethod() in {RcmdrMisc} is a wrapper i.e. inteface to deltaMethod() from {car})

  1. Effect of \(\textit{atndrte}\) at \(\textit{priGPA} = \overline{\textit{priGPA}}\):

    DeltaMethod() syntax is based on the first derivative of the regression function

\[ \textit{stndfnl}_{\,i} = \beta_0 + \beta_1 \textit{atndrte}_i + \beta_2 \textit{priGPA}_i + \beta_3 \textit{ACT}_i + \beta_4 (\textit{atndrte}_i \times \textit{priGPA}_i ) + u_i \]

with respect to \(\textit{atndrte}\):

DeltaMethod(lm2,"(b1+b4*m2)") # m2 is the average value of priGPA
##       parameter name
##     (Intercept)   b0
##         atndrte   b1
##          priGPA   b2
##             ACT   b3
##  atndrte:priGPA   b4
## 
##                 Estimate        SE     2.5 % 97.5 %
## (b1 + b4 * m2) 0.0087562 0.0026162 0.0036285 0.0139
  1. For the effect of \(\textit{priGPA}\) at \(\textit{atndrte} = \overline{\textit{atndrte}}\):
DeltaMethod(lm2,"(b2+b4*m1)")# m1 is the average value of atndrte
##       parameter name
##     (Intercept)   b0
##         atndrte   b1
##          priGPA   b2
##             ACT   b3
##  atndrte:priGPA   b4
## 
##                Estimate       SE    2.5 % 97.5 %
## (b2 + b4 * m1) 0.382031 0.078115 0.228929 0.5351

APEs/AMEs using the {margins} package

Here, we use our sample to estimate (population) the average/expected effect of changing a regressor by 1 unit (given observed covariates):

\(\textit{APE}_j = \frac{1}{n} \sum_{i=1}^{n} f(\boldsymbol{x}_i \boldsymbol{\hat{\beta}})\hat{\beta}_j\)

where \(f(.)\) is \(\frac{\partial y}{\partial x_j}\)

lm2<-lm(stndfnl ~ atndrte + priGPA + ACT + atndrte:priGPA , data=data) # repeated
summary(margins(lm2))
##   factor    AME     SE      z      p  lower  upper
##      ACT 0.0817 0.0111 7.3278 0.0000 0.0598 0.1035
##  atndrte 0.0088 0.0026 3.3469 0.0008 0.0036 0.0139
##   priGPA 0.3820 0.0781 4.8906 0.0000 0.2289 0.5351

Supervised work:

1) We expand the lm2 model to include \(\textit{priGPA}^2\) and \(\textit{ACT}^2\) as regressors:

\[\textit{stndfnl}_{\,i} = \beta_0 + \beta_1 \textit{atndrte}_i + \beta_2 \textit{priGPA}_i + \beta_3 \textit{ACT}_i + \beta_4 \textit{priGPA}_i^2 + \beta_5 \textit{ACT}_i^2 + \beta_6 (\textit{atndrte}_i \times \textit{priGPA}_i ) + u_i\]

Estimate the model
# Fill-in and uncomment the following lines
# lm4 <-lm(stndfnl ~ ... , data=data) 
# summary(lm4)
Test \(H_0: \beta_1 = \beta_6 = 0\)
# Fill-in and uncomment the following lines
# ?linearHypothesis
# linearHypothesis(lm4, ... )

For results verification, see Wooldridge, Introductory econometrics, Example 6.3

2) Use R (DeltaMethod) to calculate marginal effects and standard errors. Use mean values of all regressors as your “interesting” observations.

# 
m3 <- mean(data$ACT)
# 
# Effect and s.e. of atndrte
#
# Effect and s.e. of priGPA
#
# Effect and s.e. of ACT
#

3) Use R (margins) to calculate APEs of all regressors, interpret

# APEs 
#

For additional examples, see margins: https://cran.r-project.org/web/packages/margins/vignettes/Introduction.html.